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Abstract. Tiie problem of tiie determination of tlie cliarge density from limited in- 
formation about the charge form factor is an ill-posed inverse problem. A Bayesian 
probabilistic approach to this problem which permits to take into account both errors 
and prior information about the solution is presented. We will show that many classical 
methods can be considered as special cases of the proposed approach. We address also 
the problem of the basis function choice for the discretization and the uncertainty of 
the solution. Some numerical results for an analytical model are presented to show the 
performance of the proposed method. 

1. Introduction 

Elastic electron scattering provides a mean of determining the charge density of a 
nucleus, p{r), from the experimentally determined charge form factor, F{q). The 
connection between the charge density and the cross section is well understood 
and in plane wave Born approximation F{q) is just the Fourier transform of p{r) 
which for the case of even-even nuclei, which we shall consider, is simply given by 

/•OO 

F{q)=47T Jo{qr)p{r)dr (1) 

JQ 

where Jq is the spherical Bessel function of zero order and q is the absolute value 
of the three momentum transfer. Given that the experimental measurements are 
performed over a limited range at a finite number of values of the momentum 
transfer q, a unique determination of p{r) is not possible since the resulting inverse 
problem is ill posed. 

One of the generally accepted procedures for determining p{r) is to expand it 
in a basis and then determine the expansion coefficients from a least squares (LS) 
fit to the experimentally measured values of F{q) |jl|, ||, |[ [l^. The following 
questions then arise: how to choose a basis and how to determine the order of the 
expansion? Another problem with the LS methods is that increasing the number 
of terms in the expansion generally leads to non physical oscillations in the charge 
density in spite of the fact that the charge form factor is well reproduced at the 
experimentally determined values oi q |@, |l^. Finally, due to the fact that the 
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problem is inherently ill posed, a small error in the data (experimental errors 
or measurement noise) will produce large variations in the solution which is not 
acceptable in practical situations. 

What we are going to do is to show how a Bayesian approach can be helpful to 
give both correct and reasonable answers to the aforementioned questions and to 
propose new methods which are more stable with respect to the errors and finally 
to give procedures to put the correct error bars on the proposed solutions. 

2. Fundamentals of the Bayesian approach 

Let us start by discretizing the problem in the usual manner by expanding p{r) in 
a basis 4'n{'r)'- 

^ \ r>Rc ^ ' 

and substituting it in yields 

F{q) = 47r / Jo(g?-) V W dr 

= 47r^a„ / Mqr)(j>n(r)dr (3) 

n=l Jo 

Now, defining 

rR. 

^m,n = 47r / Jo(g„r)0„(r)dr (4) 







we obtain 

Aa + e (5) 

where a is a vector containing the coefficients {a„, n = 1, • • • , iV}, Fc is a vector 
containing the form factor data {Fc{qm), m — 1, • • • , M} and A an (M x N) matrix 
containing the coefficients Am.n given by (^. Note also that when the vector a is 
determined, we can calculate p = {p{rk), fc = 1, • • • , K} by 

p = *a (6) 

where ^ is a. K x N matrix with the elements $fe„ — (pnirk)- 

The vector e is added to take account of the errors in both measurement noise 
and due to discretization. We assume that the components {cm,™ = 1, • ' ' t 
are additive, zero mean (no systematic error), mutually independent (no correla- 
tion) and independent of a, and they can only be characterized by their common 
variance (7^ . This hypothesis is reasonable unless we know more about its charac- 
teristics. 

Note that we have not yet discussed the choice of the basis functions (/>„ and 
the determination of the expansion order, N. We will come back to these questions 
later. Let us now see how the Bayesian estimation approach works. 
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The main idea behind the Bayesian probabihstic approach is to represent the 
uncertainty or any lack of knowledge or any diffuse prior knowledge about a quan- 
tity by a probability law. For example, the knowledge (or the hypothesis) that 
{£rm m = 1, • • • , Af} are zero mean, mutually independent and that they are only 
characterized by their common variance can be described by choosing a Gaus- 
sian probability distribution for them. One may also use the Maximum Entropy 
Principle to enforce this choice. This means that we can write 



or p{e„i) = A/'(0, o-^ ) and 



m— 1 ^ ' V. € 

or simply p{e) = Af{0, cril) where I is the (A/ x M) unitary matrix and is the 
common variance of e,„ for all m. Now, using this model (^), we can define the 
conditional probability law 

p(Fc|a)=(^^^cxp{-^||F.-Aa|p} (9) 

or p{Fc\a) = N{Aa^ '^ll)- It is usual to call p{Fc\a) or its logarithm, considered 
as a function of a, the Likelihood. 

One can stop here and define the solution of the problem (||) as the vector a, 
which maximizes the likelihood (ML): 



2 = argnmx{p(Fc|a)} (10) 



or equivalently 



a = argmin {— lnp(i^c|a')} (H) 
which, in the case of a Gaussian distribution (^) becomes 

a = argnfin{||Fc - Aa|p)} (12) 

and we find here the LS solutions given by: 

{A^A)a^A^Fc (13) 

The main problem with these solutions is that, very often, the matrix A^ A is 
either singular or at least ill-conditioned, so that the solutions are very sensitive 
to the errors in the data or even on the round-off errors during the numerical 
calculation. 

The Bayesian approach can do better. In fact, before looking at the data, we 
may have some prior knowledge about a. For example, we know that a € and 
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that we might prefer those vectors which have a minimal norm ||a||. To translate 
this prior knowledge and this preference, we may assign a Gaussian probability 
distribution to the vector a: 

P(«)-(^^^-p{-^ll«ll^} (14) 

where gives an idea about the scale of the norm of the vector a. Now, using 
the Bayes rule we can calculate the posterior probability law 

, , , p(Fc\a)p(a) , , 

pi-\F.) = (15) 

where the denominator 

= J p{F,\a)pia)da (16) 

is the normalization constant (called sometimes the Evidence) [^ . 

This posterior law contains all the information we may wish about the solution. 
For example, we may want to know what is the probability that a < a <a. This 
can be calculated by 



P{a<a<a)= j p{a\Fc)da 
J a 



(17) 



Or, we may be interested only in one of these parameters a„ and want to know 
what is the probability that a„ < a„ < a„. This can be calculated by 



P(a„ < a„ < a„) = / p(a„|Fc)da„ (18) 

where the marginal posterior law p(a„|i^c) can be calculated by 

p(a„|Fc) = y" • • • y"p(a|Fc) dai • • • da„_i da„+i • • • dflAT. (19) 

We can also simply define as the solution the vector a which corresponds to the 
mean value of the posterior law -called Posterior mean (PM) estimator: 

a = j ap{a\Fc)Aa (20) 

or the vector a which maximizes this posterior distribution -called Maximum a 
posteriori (MAP) estimator: 

a — argmax{p(a|Fc)} ~ argmin {— lnp(a|i^c)} (21) 



4 



or even the vector a whose components a„ correspond to the maximizer of the 
marginal posterior law (^9|) -called Marginal MAP estimator: 

a„ = argmax{p(a„|Fc)} = argmin {- lnp(a„|Fc)} . (22) 

In the following wc consider only the MAP estimator (pl|). Using the proba- 
bility distributions (||) and ( p^ ) in (p^, the MAP solution is given by: 

a = argnnn{||Fc - a|p + A||a|p)} (23) 

where A = [fr^jaaf' and we find the minimum norm least squares (MNLS) solution 
which is given explicitly by 

a = (A*A + A7)-iA*Fc. (24) 

Comparing ( |l^ ) and ( p^ we see that, for a given N ^ the matrix (A*A + A/) is 
always better conditioned than the matrix (A*A) and so the solution (|2^ ) always 
is more stable than the solution (p^). 

We may also want some information about the uncertainty of this solution. For 
this we can use the posterior law p(a|Fc)- For example, using the likelihood (|) 
and the prior law (p^, it is easy to show that the posterior law is Gaussian, i.e. 
p(a|Fc) = A/'(a,P) with a given by (||) and P = {A^ A AJ)-^ We can then 
use the diagonal elements of the posterior covariance matrix P to calculate the 
posterior variances of the estimates, i.e. Var(a„) = P,m and so put the error bars 
on the solution. When the posterior law is not Gaussian, we can always calculate 



E(a„) = y a„p(a„|Fc)da„ (25) 

and 

Var(a„) = ^ (a„ - E(a„))^ p(a„|Fc) da„ (26) 

but in general we may not have explicit expressions for these integrals. We can 
however do numerical calculation either by approximating the posterior law by a 
Gauusian law or by a stochastic integral calculation. 

One question still remains: how to determine A and iV? Three approaches are 
possible: 

1. Assign them experimentally from the data using some knowledge on the 
physics of the problem. For example, the Parseval-type relation between p{r) 
and Fc(g): 

y 4^r2p2(^) dr = ^ y 4^q2^c(q) (27) 
can be used to estimate a\ by: 

N M 



n— 1 m— 1 

and having an estimate of the noise variance we can determine A. 
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2. Consider A and N as two extra parameters (hyper parameters) to estimate 
jointly with the unknown parameter a. We can then assign a prior law for 
them. For example Jeffrey's priors p(A) ~ j for A and a uniform p{N) = 
1 /N^^^ for N. (Other choices are possible, for example a Gamma prior A which 
eliminates A = and X = oo and a binomial prior for TV which eliminates 
N = andN = 7V„,,.) 

Finally, we can estimate them jointly with a by 

(a,A,iV)=arg max {p(a, X, NlFc)} (29) 
(a,x,N) 

where 

p{a,X,N\F,)^piF,\a,X,N)p{a\X,N)p{X)p{N). (30) 

We must however be careful to verify that this joint criterion has at least a 
local optimum. 

3. Consider A and TV as two extra parameters as in the precedent case but not 
on the same level. This means that we can try to estimate them first by 

(A, TV) = arg max {p(A, TV|F,)} (31) 

(A, AT) 

where 

p{X,N\Fa) = J p{a,X,N\Fc)da (32) 

and then use them in (|T]). Note, however that finding an analytical expression 
for p(X, TVji^c) is not always possible and its numerical calculation may be very 
costly. 

4. Consider A and TV as two nuisance parameters, integrate them out and esti- 
mate a directly by 



a = argnmx{p(a|Fc)} = argniax 1^ J p{a, X,n\Fc) dX^ 



(33) 



(For more details on these methods, their relative characteristics, their practical 
implementations and their relatives performances see p^). 
Still one question remains: the choice of the basis-functions. 



3. Choice of the basis-functions 

Two approaches are used to select the basis functions. We call them the operator 
based parametric approach and the non parametric approach and we will discuss 
both in detail in the following sections. We propose then a new third approach 
which tries to eliminate the limitations and to keep the advantages of the previous 
approaches. We call this third approach physically based parametric. 
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3.1. Operator based parametric approach 

The first approach is to choose special purpose basis functions based on the 
properties of the operator linking the data to the unknowns. For example in our 
case, due to the fact that the kernel of the integral operator of the direct problem 
is a Bessel function, we may also use the Bessel functions as the basis functions 
for p{r) 

= / E^=i ank{qnr) r < 
^ \ r > Rc 

This will permit us, using the orthogonality relation 



L 



Ji{qnr) Ji{qrar) dr ^ Ji+iiQnRc) Sn,7n, (35) 



to find an explicit expression for the charge form factor as a function of the coef- 
ficients a„ : 

F{q) = ^ «" t u 2 f ^2 sin 9^ • 36 



With this choice, note also that, if the form factor F{q) was known exactly at 

n-TT 



qn = jf- then the coefficients a„ could be calculated analytically by 



Fjqn) ^o7A 

a„ = ^. (37) 

27ri?3[Ji(q„i?)]' 

In general, however, the cross section is measured at momentum transfers dif- 
ferent from q-n — 

Now, assume that we are given M measurements at arbitrary momentum trans- 
fers q — {qi, q2, ■ ■ ■ , qn} and we wish to determine the N expansion coefficients 
a — {oi, 02, . . . , aj^}. In this case Eq. ( ^6| ) leads to 

Fc = Aa + e (38) 

as in (||). 

The main advantage of this approach is the fact that a is a small dimension 
vector and so is the matrix A and we have an explicit analytical expression for 
calculating its elements. 

But at least one main disadvantage to such a choice is that our prior knowledge 
on a may be limited. For example, if we know that p(r) is a positive function we 
cannot easily incorporate this information in the parameters, a. 

3.2. NON PARAMETRIC APPROACH 

The second approach is to choose the basis-functions as general as possible and 
independently of the direct problem operator, for example, either: 

Mr)^Sir~nA) (39) 
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or 

1 if (n - 1)A < r < nA 



'^"^'^) = \ elsewhere ^^^^ 

with A chosen appropriately small (maximum needed resolution) to be certain 
we are able to approximate any function p(r) as precisely as desired. But, this 
means that TV will probably be large. This is a disadvantage, but this can be 
compensated, as we will see further, by the fact that the coefficients a„ now have 
a direct physical meaning: the samples of p(r) in the first case and the mean values 
of p{r) in the intervals (n — 1)A < r < nA in the second case. This means, for 
example, that the prior knowledge such as the smoothness or the positivity of the 
function p{r) can be transmitted to the coefficients a„ easily. 

Let us choose ( p^ and go further into the details. Replacing (|^) with the 
basis-functions (391) in dll) we obtain: 



F{q)^^an drinr'^ Jo{qr)S{r ~ nA) ^ 4:n^anA{nAf JoinAq) (41) 

n=l n=l 

Denoting by 

Am,n = 47r A {nAf Jo(nAg„0 (42) 



we obtain Fc = Aa + e as in (|38|). If we use (^) in place of (pOD, the only change 
will be in the expression of the Matrix elements Am,n which become 

/•nA 

^m,« = 47r/ Jo{q„,r)<:\r. (43) 

J(ri-1)A 

To make a distinction between this approach and the preceding one, let us 
denote a by p and Ahy B: 

F, = Bp + e (44) 

Let us now compare ( p^ ) and (p8|): a in (^) is a vector of small dimension 
while p in ( |4^ ) is a vector of much larger dimension. 
Here we can also define either the LS solution: 

p = argmm{||F,-Bpf} = {B'B)-^B'F^ (45) 

or the MNLS solution: 

p = argmm{|lF,-Bp||2 + A||p||2} = {B* B + XI)-' B* F^ (46) 

but neither of these solutions may be satisfactory. 

In ( ^ ) it is possible to incorporate the smoothness and the positivity of the 
function p(r) into a more appropriate prior distribution for the components p„. 
For example, to enforce the smoothness of p{r) we can assign 

N 

P{P) = P{Pl)Y[pipn\Pn-l) = J^{pn-l,(Tl) (47) 
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with 



p{pi) =7V(/9o,CTo) cx exp 



p{Pn\p) =p(p„|p„-i) =7V(p„_i,cro) CX exp 
which leads to 



p{p) oc exp 



1 ^ 



Using this prior distribution in ([l5|), the MAP estimator becomes 



with 



p — argniax{p(p|i^c)} — argniinj J(p)} 



N 

J{p) = \\F,-Bp\\^+Xj2iPn-Pn-l) 
n=l 



(48) 
(49) 

(50) 
(51) 
(52) 



Defining the matrix 



D 



(\ -1 
1 



V 



(53) 



it is easy to show that 

J(p) = ||F,-Bp||2 



1 -1/ 



A||(r»p||2 + A(pi-po)' 



(54) 



Let us temporarily assume that p\ — po- We then have an explicit solution for the 
minimizer of (E3) which is given by: 



p = {B*B + \D^D)-^B^Fc 



(55) 



Comparing this solution with the MNLS solution ( pq ) gives us the possibility to 
see the difference in which the term D^D is used in place of /. Indeed, due 
to the fact that D corresponds to a first order derivative, we may designate the 
MNLS solution as the zero order regularized solution in contrast to the first order 
regularized solution. It is possible to extend this to more general regularized 
solutions by an appropriate choice of the matrix D. 

Now, let us go back to (p^). pQ is now a new extra hyper-parameter which 
may play a great role in the solution of our inverse problem where the data do not 
contain information about the DC level of the function p{r). 

One way to enforce the positivity of the solution is to choose a prior distribution 
such as : 



P{P) 



N 

E 

n=l 



P{Pn 



N 



^r(a) 

n— 1 ^ ^ 



exp[-/3p„ 



(56) 
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which can also be written as: 



p(p) oc exp 



N 



^(1 - a) Inpn + f3pn 



(57) 



and which is called an Entropic prior in |2C[| . Using this prior law in (|15|), the 
MAP estimator becomes 



p = argmax{p(p|Fc)} = argniin{ J(p)} 



(58) 



with 



N N 

J{p) = - Bpf + Ai ^ lnp„ + A2 ^ p„ (59) 

n—1 n—1 

where Ai and A2 are related to a, /3 and a^. Other choices are possible pOt - 
To enforce both positivity and the smoothness we propose here to choose 



p{p) (x exp 
which leads to 



N N N 

-Ai ^{Pn - Pn-lf - A2 ^ \npn ^ A3 ^ p„ (60) 
n—1 n—1 n—1 



p = argmin |J(p)} (61) 



with 



N 

J{p) = - Bp\\^ + J2 MPn - Pn-lf + A2 lnp„ + AaPn (62) 

n=l 

3.3. Physically based parametric approach 

In this approach we choose special purpose basis functions based on the physics 
of the problem. For example, in our case, since the charge density is a single- valued 
function defined in a finite domain, the Fourier-Bessel (FB) basis functions which 
satisfy both the orthogonality and the concentration property conditions, can be 
used for expansion: 



p{r) 







r>Rr 



(63) 



where Qn = Excepted the original motivation, this choice is exactly the same 
as in the first approach and all the relations developed and discussed there can be 
used. 

We may choose other basis functions which are more appropriate to translate 
our prior knowledge on the desired solution. For example, in our case, we know a 
priori that, the solution is smooth, positive and a decreasing function. Then we 
can choose the following function : 



p{r) 



Y.n=i exp(-g„r2) r < Rc 
r>Rr 



(64) 
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Using this expansion in cq. (|^) we find Fc = Aa + e as in (||) or in , where 
Am,n = 4:Tr drr^ Jo((?r„?')exp(-g„r^). (65) 

With this choice we keep the main advantage of the first approach which is the 
small dimension of the vector a and the main advantages of the second approach 
which is the translation of our prior knowledge of the positivity of the function 
p{r). This is due to the fact that if we impose the positivity constraint on the 
coefficients a„ we insure that the solution remains always positive. 

In the next section we will illustrate the performance of these different solutions 
for the estimation of the charge density from elastic electron scattering data. 

4. Issues on the uncertainty of the solution 

In any scientific problem solving, a proposed solution should be given in any way 
with a measure of its uncertainty or confidence. In Bayesian approach, the poste- 
rior probability gives us naturally the necessary tool. To see this, let come back 
to our problem and make a summary. We have a set of data F and we want 
to estimate p{r) or more precisely p for some locations r^. Let assume that we 
have chosen a constant discretization step and so, we want to estimate a vector 
P = {p{n)i = 1,---,K}. 

We presented two approaches: parametric and non-parametric. In the first 
approach, we have 

F = Aa + e (66) 
p = *a (67) 

and in the second 

F = Bp + e (68) 

In both cases, we are interested to p. In the first approach, we assigned p{a) and 
p{F\a), calculated the posterior p(a|i^), defined a solution a for the parameters a, 
and finally, a solution p = $2 for p. In the second, we assigned directly p{p) and 
p{F\p), calculated the posterior p(p|F), and finally defined a solution p. In both 
cases, we can use the posterior laws to quantify the uncertainty of the solutions. 
There are, at least, three approaches: 

— Simply generate samples from the posterior law p(p\F) using for example a 
monte carlo method, and show all these samples to see the distribution of the 
proposed solution. 

— Calculate the posterior mean and the posterior variance of the solution at each 
point either analytically (when possible) or numerically using for example the 
samples generated by a monte carlo method. 

— Calculate the posterior mean and covariance of the solution either analytically 
(when possible) or approximate it numerically by any quadrature algorithm. 
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To illustrate this, let eonsider the cases where all the probability laws are 
Gaussian. Then, all the calculations can be done analytically. The following 
summarizes all the steps for the calculation of the solutions and their posterior 
covariances in the above-mentioned two cases: 



Non-parametric 


Parametric 


F = Bp + e 


F = Aa + e 
p = *a 


p{p) = N{po,crlPo) 


p(a) = M{ao,all) 


p{F\p) = N{Bp,all) 


p{F\a) = N{Aa,aU) 


p{p\F) = N{p,Pp) 


p{a\F) = J\f{a,Pa) 


p = [B'B + XP-Y'B\F-Bp,) 


a = [A^A + \I]-^A'{F - Aao) 


Pp = [B'B + XPoT' 


Pa = [A'A + XP^T^ 


with A = CTe /ffp 


with A = erf 






p{p\F) = AA(*S,*Pa**) 




Special cases 




Special cases 




"p= [B*B]-^B*F 
Pp = 


a, ^0 


.Pp-0 


ap^O < 


'P = Po 
.Pp = 


o-a — >0 


p = *ao 
.^^^ = 


CTj > OO < 


'P = Po 

Pp = ^iPo 


CTc > OO < 


p = *ao 


ap — > OO < 


'p = [B^B]^B*F 
Pp = al[B'B]-^ 


CTa > OO < 





Table 1: A comparison between parametric and non-parametric approaches. 
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When the posterior covariance matrix Pp is calculated, we can use it to give 
some information about the uncertainty of the solution. For example, we can 
use its diagonal elements to calculate at = y/Pkk and use it to error bars on the 
solution. 



5. Numerical experiments 

In order to demonstrate the preceding considerations we make use of the following 
analytical model. For a charge density given by a symmetric Fermi distribution g 

cosh(i?/d) 



(69) 



cosh(i?/(i) + cosh(r/d) 

an analytical expression for the corresponding charge form factor can easily be 
obtained §, §: 

47r^ ad cosh (i?/c?) \ R cos{qR) nd sm{q R) cosh{Tr qd) 



F{q) 



q smh{R/d) 



smh{TTqd) 



sinh (nqd) 



(70) 



Only two of the parameters a, R and d are independent since the charge density 
must fulfill the normalization condition 



47r J p[r) dr = Z. 



(71) 



Figure 1 shows the theoretical charge density p{r) of ^^C (Z=6) obtained from 
(H) for r € [0, 0.8] with R = 1.1 A* and d = 0.626 fm and the theoretical charge 
form factor Fc{q) obtained by ( [70| ) for q G [0,8] fm^-^ and the nine simulated 
experimental data: 

q = [0.001, .5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0] fm"^ 




Figure 1. Theoretical charge density p(r) [left], charge form factor log ||Fc(q)|j and the data 
[stars] used for numerical experiments [right]. 
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5.1. Experiments with operator based parametric models 

We use these data in the parametric model with i?c = 8 fm and estimate 
the coefficients a by 

LS : a = {A^A)-^A^Fc 
MNLS : {A^A + Xiy^A^Fc 

MAP : a = [A^ A + XD^Dy^A^F^ 

Then using these coefficients we calculate p{r) by (|3^ ) and Fdq) by (^6|). 

Figure 2 and Figure 3 show the reconstructed charge densities p and the corre- 
sponding charge form factors Fc obtained by LS and by MNLS for iV = 5 and 
N = 10. 

Figure 4 shows the reconstructed charge densities by LS and by MNLS for different 
expansion order N from 5 to 10. Note that the LS solutions are very sensitive and 
vary greatly with N , but the MNLS solution stays more stable with respect to N . 

Figure 5 shows the reconstructed charge densities and the corresponding charge 
form factors obtained by MNLS and MAP for N = 30. 




Radius Momentum transfer 



Figure 2. Parametric reconstruction of p(r) obtained by LS (point) and by MNLS (dotted) for 
N = 5 [left] and the corresponding reconstructed charge form factors log 11^0(9)11 [right]. Two 
solutions are practically indistinguishable and both not very satisfactory due to a large bias of 
the solution for small radius r. Note also that both solutions fit well the data. 
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Charge density 4 




01 234567801 2345678 

Radius Momentumtransfer 



Figure 3. Parametric reconstruction of p{r) obtained by LS (point) and by MNLS (dotted) for 
A'' = 10 [left] and the corresponding reconstructed charge form factors log ||_Fc(g)|| [right]. Note 
that the LS solution fits very well the data but is very unstable but the MNLS solution, which 
does not fit perfectly the data, is at least more stable. 

Charge density Charge density 

3| i i i i i i i 1 1.2| i i i i i i i 

2 -■ ■ ■ - 1 - '■ 




Radius Radius 



Figure 4. Parametric reconstruction of p(r) obtained by LS (left) and by MNLS (right) for 
different values of W := 5 : 1 : 10. Note that the LS solutions vary greatly with A'', but the 
MNLS solutions stay more stable with respect to TV. 
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012345678012345678 
Radius Momentum transler 



Figure 5. Parametric reconstruction of p(r) obtained by MNLS (point) and by MAP (dotted) for 
N = 30 [left] and the corresponding reconstructed charge form factors log ||i^c('j)|| [right]. In this 
case the LS solution is completely unrealistic and is not presented. The MNLS solution has a 
large bias for small raduis. The MAP solution is very satisfactory. Note also that both solutions 
satisfy the data constraint practically in the same way. 



5.2, Experiments with non parametric models 

The same data are then used with the non parametric model (|4^ ) with N = 100 
and i?c = 8 fm and p is calculated by 



LS : 


P = 


(B'B) 


-^B'F, 


MNLS : 


P = 


(B'B 


^XI)-^B*Fc 


MAPI : 


P = 


(B'B 


^ XD^D)-^B*Fc 



Figure 6 shows the estimated p and the corresponding Fc by MNLS and by 
MAP. 

Figure 7 shows two solutions obtained by a parametric and a non-parametric 
method and their associated error bars. 
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Figure 6. Non parametric reconstruction of p{r) obtained by MNLS (point) and MAPI (dotted) 
for N = 100 [left] and the corresponding reconstructed charge form factors log 11^0(9)11 [right]. 
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Figure 8. Uncertainty in parametric and non parametric methods: 
Left: Parametric reconstruction of p(r) obtained by MAPI 
Right: Non-parametric reconstruction of p(r) obtained by MAPI 
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6. Conclusion 

We considered the problem of the determination of the charge density from a 
limited number of charge form factor measures as an ill-posed inverse problem. We 
proposed a Bayesian probabilistic approach to this problem and showed how many 
classical methods can be considered as special cases of the proposed approach. We 
addressed also the problem of the basis function choice for the discretization and 
the uncertainty of the solution. We illustrated the performances of the proposed 
methods by some numerical results. 
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